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by 
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The Boeing Company, Seattle, Washington 


1. 0 ABSTRACT 

Analytical methods based on linear theory are presented 
for predicting the thermal stresses in and the buckling of 
heated structures with arbitrary uniform cross section. The 
structure is idealized as an assemblage of laminated plate-strip 
elements, curved and planar, and beam elements. Uniaxially 
stiffened plates and shells of arbitrary cross section are 
typical examples. In such structures, the common practice is 
either to "smear" the 3tiffeners or to model them as beams along 
their elastic axis. In contrast, the idealization used here 
retains the relative stiffnesses and restraints between parts 
of each stiffener. This may have significant effect on the 
distribution of thermal stresses and the buckling 
characteristics. For the buckling analysis the structure or 
selected elements may be subjected to mechanical loads, in 
addition to thermal loads, in any desired combination of inplane 
transverse load and axial compression load. The analysis is 
also applicable to stiffened structures under inplane loads 
varying through the cross section, as in stiffened shells under 
bending. The buckling analysis is general and covers all modes 
of instability. The analysis has been applied to a limited 
number of problems and the results are presented. These while 
showing the validity and the applicability of the method do 
not reflect its full capability. 



2 . 0 SYMBOLS 
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B. , 
ID 
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ID 


length of the structure 

temperature distribution coefficients, 
equation (B.3) 

extensional stiffnesses i, j = 1, 2, 6 
equation (A. 7) 

cross-sectional area of beam element 

developed width of plate-strip element 

stiffnesses (i, j = 1, 2, 6) associated 
with coupling between bending and extension, 
equation (A. 8) 

diagonal matrix of longitudinal variables, 
equations (B.13) , (C.4) and (C.5) 

transfer matrix, equations (C.ll) and (C.12) 

displacement vector of elements 

displacement vector of sub-elements 

displacement vector of total structure, 
equations (B.39) and (C.20) 

bending stiffnesses (i, j = 1, 2, 6) equation 
(A. 9) 


E ll r E 22 Young's modulii of orthotropic material 

f force vector of elements 



f 



F 



equivalent thermal force vector of plate- 
strip elements, equations (B.23), (B.28) 
and (B.21) respectively 

force vector of sub-elements 

equivalent thermal force vector of beam 
elements, equation (B.34) 

equivalent thermal force vector of total 
structure, equation (B.39) 

shear modulii of orthotropic material 
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distance to the k t " lamina from the reference 
plane 

polar moment of inertia of beam element 
moments of inertia of beam element 
torsion constant for beam element 
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effective torsion constant for heated beams, 
equation (A. 56) 

diagonal matrices spring constants 
spring constants 

coefficients of characteristic equations 
(B. 10) 

particular solution coefficients, equation 
(B.7) 
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M ll' M 12' M 22 
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number of laminas 

displacement ratio coefficients, equations 
(B.ll) 

linear differential operators, equations 
(A. 18 to A* 23) 

Fourier coefficient number for thermal 
stresses and axial half-wave number for 
buckling, equations (B.2) and (C.3), 
respectively 

moment resultants, equation (A, 6) 

torque on the beam element 

moment resultants in heated beam elements, 
equations (A, 54) and (A. 55) 

moment resultants in heated plate-strip 
elements, equation (A, 11) 

moment resultant, equation (A. 25) 

stress resultants, equation (A, 5) 

force resultants in heated plate-strip 
elements, equation (A, 10) 



total applied inplane loadings 


applied mechanical inplane loads, equation 
(C.l) and (C.2) 

applied inplane preloads, equation (C.l) 
and (C.2) 

effective stress resultant in y-direction, 
equation (A. 26) 

displacement parameters, equation (B.6) 

buckling displacement parameters for sub- 
elements, equation (C.3) 

axial load in beam element induced by axial 
displacements equation (A. 52) 

buckling load of plate-strip elements, when 
the longitudinal (x) sides are completely 
restrained 

axial load in heated beam element, equation 
(A. 53) 

applied axial load in beam 

lateral shears on the beam element equation 
(A. 49) and (A. 51) 

orthotripic material constants (i, j = 1, 

2, 6) , equation (A. 2) 

effective transferse shear parallel to z- 
axis in plate-strip elements, equations 
(A. 24) and (A. 36) 

reference plane radius of the curved plate- 
strip element 

elements of coefficient matrix R, equation 
(B. 9) 

stiffness matrices for plate-strip elements 


stiffness matrix for beam element 



reduced stiffness matrices for plate-strip 
element, equations (B.23), (B.28), (C.17) 
and (C.18) 

stiffness matrix of total structure 

thickness of layer of a laminate 

temperature distributions (from ambient 
conditions) 

transformation matrices 

effective inplane shear in plate-strip 
elements, equation (A. 27) 

displacements at the reference plane of 
plate-strip element and at the shear center 
of beam element 

buckling displacements at the reference 
plane of sub-elements 

displacement coefficient for beam element, 
equation (B. 33) 

displacement coefficients for plate-strip 
element, equation (B.6) 

buckling displacement coefficients for sub- 
elements, equation (C.3) 

orthogonal coordinates, figures (A. 2) and 
(A. 5) 

matrices for plate-strip elements, equations 
(B. 13) to (B.16) 

matrix for beam element defined in equation 
(B.34) 

matrices for sub-element, equations (C.4) , 
(C.5) 

plate-strip element offsets, Figure 25 

distances measured parallel to the y and 
z axes respectively, from the shear center 
of the beam element to its neutral axis 

displacement parameter, equations (B.6) 
and (C.3) 
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Subscripts 

AD 


b 

BC 


(g) 

G 

i 

k 


coefficients of thermal expansion 

displacement parameter, equation (B*6) 

buckling displacement parameter, equation 
(C.3) 

shearing strain 

warping constant of beam element 
normal strains 

angle subtended at its center of curvature 
by the plate-strip element. Figure 25 

rotation or twist 

beam twist coefficient, equation (B.33) 

change in curvatures 

Poisson's ratio 

normal stresses 

shear stress 

transformation angles, equations (A. 40) 
to (A. 42) 

angle between y G axis and the chord of 
plate-strip element. Figure 25 

quantities along the side AD of plate-strip 
element 

quantities related beam element 

quantities along the side BC of plate-strip 
element 

index for sub-element number 
quantities related to global axes 
index corresponding to characteristic roots 
layer index 



p 


particular solution parts 
quantities with offset effects 


s 

A subscript preceded by a comma indicates partial differentiation 
with respect to the subscript. 

Superscripts 

k layer index 

T matrix transpose 

° quantities in the reference plane of plate- 

strip element 

+ quantities along 

strip element 

quantities along 
strip element 


the side y = + - of plate- 


the side y = - - of plate- 



3.0 INTRODUCTION 


Aerospace vehicles and missiles often operate under severe 
thermal conditions, thus subjecting the structural components 
to thermal gradients. The resulting thermal stresses can have 
significant effects on the characteristics of these structures. 
Refs. 1, 2 and 3. 

Being a self-equilibrating system, the thermal stresses 
are both tensile and compressive. In thin-walled structures 
like stiffened plates and shells these stresses reduce the 
structural performance under mechanical loads and precipitate 
buckling, Refs. 4 and 5. Thermal stresses can also cause 
buckling by themselves, if the temperatures are sufficiently 
high. The usual approach to buckling analysis treats the thermal 
stresses as pre-loading in conjunction with the established 
methods of isothermal buckling analyses. Refs. 1 to 6 are 
typical examples. Thermal stress analysis is thus an essential 
first step to the proper understanding of the buckling behavior 
of these structures under combined thermal and mechanical loads. 

The analysis of Ref. 7 for buckling under mechanical loads 
is extended here for thermal stresses and buckling of laminated 
composite structures of uniform cross section like stiffened 
panels subjected to these stresses and mechanical loads. As 
in Ref. 7 the structure is idealized as an assemblage of 
laminated plate-strip elements, curved and planar, and beam 
elements. Each element extends the full length of the structure. 
Element edges normal to the longitudinal axis are assumed to 
be simply supported with no restriction on the axial (warping) 
displacement. 

The analysis of Ref. 5 is primarily intended for complex 
shells of revolution. It can also be applied to stiffened 
panels, retaining the discreteness of the stiffeners, by 
considering these panels as ring stiffened shells of very large 
(infinite) radius. Each ring stiffener is treated as a single 
equivalent beam. Such an idealization is justified when 
considering the general behavior of complex shells of revolution. 
However, when applied to stiffened panels the idealization of 
stiffeners in Ref. 5 is insufficient to predict any local modes 
of buckling in the stiffeners or any coupled modes involving 
a contribution from local deformations in the stiffener. 

The present analysis enables each stiffener to be idealized 
as an assemblage of the types of elements mentioned before. 
Intuitively, this is a better representation of the true behavior 
of the stiffener. Further, it is possible to distinguish between 
riveted and bonded connections by means of the modelling 
techniques discussed in Refs. 7 and 8. 



The curved and flat plate-strip elements are in general 
laminated. Linear or parabolic thermal gradients may be 
specified across the width of individual elements. An arbitrary 
thermal gradient through the thickness is approximated by a 
series of layers, each at a constant temperature. The 
temperature is constant in the axial direction. The stress- 
strain equations for each lamina assume orthotropy with respect 
to the longitudinal and the transverse axes of the element. 

Elementary theory of bending and torsion of heated beams 
under axial loads is used for the beam elements. The section 
properties of laminated beams are calculated in an approximate 
manner. Laminated circular and rectangular beam elements are 
considered in detail. For these, the temperature is arbitrarily 
constant in each lamina and in the axial direction. For the 
latter type of beam element, either identically linear or 
parabolic thermal gradients are allowed across the width of 
each lamina. 

The analysis considers offsets between elements and the 
effects of arbitrary elastic restraints along any external 
longitudinal side (i.e., not connected to other elements) of 
the flat or curved plate-strip elements. 

In the thermal stress analysis, a truncated Fourier sine 
series is used to approximate the constant temperature in the 
longitudinal direction of each element. Thermal stresses 
evaluated for each harmonic of the Fourier series are superposed 
to obtain the resultant thermal stresses corresponding to the 
boundary conditions discussed before. 

The buckling analysis treats the axial thermal load in 
each beam element as preloads. The above thermal stress 
resultants, while varying across the width of each plate-strip 
element, are assumed constant in the longitudinal direction. 
Because of this transverse variation of the stress resultants, 
the governing differential equations for each plate-strip element 
have variable coefficients. Thus, each such element is sub- 
divided and consequently the variable-coefficient differential 
equations are reduced, Ref. 9, to corresponding sets of 
equivalent constant-coefficient differential equations. 

The mechanical loading is basically longitudinal 
compression. However, as in Ref. 7, selected plate-strip 
elements may also be loaded in the inplane transverse direction. 

A typical example is biaxially loaded stiffened panels under 
thermal gradients. In such panels, it is reasonable to assume 
that the inplane mechanical loading transverse to the stiffener 
direction is carried entirely by the plate-strip elements in 
the plane of the skin. 



The analyses of thermal stresses and buckling satisfy 
compatibility on the inter-element boundaries by considering 
four kinematic degrees of freedom. These correspond to the 
three translational displacement vectors and the rotational 
displacement vector about the longitudinal axis. Element 
stiffness matrices relating the above displacement vectors to 
the corresponding force vectors on the inter-element boundaries 
are appropriately merged to yield a set of non-homogenous 
equations for the thermal stress analysis. For the buckling 
analysis a similar procedure gives a set of homogeneous 
equations. The buckling load is obtained by iteratively solving 
the resulting non-linear eigenvalue problem in a manner analogous 
to Ref. 7. The eigenvector solution determines the buckling 
mode shape. The buckling analysis is general and as shown in 
Ref. 7 yields the lowest buckling load irrespective of the type 
of buckling. 

The basic assumptions in the analysis are: 

(a) The material is linearly elastic. However, the elastic 
properties corresponding to the average temperature 

in each lamina may be used. 

(b) Each lamina is orthotropic. 

(c) The Kirchhoff-Love hypothesis is used for the 
deformation across the thickness. 

(d) Thermo-elastic coupling is ignored. 

(e) Pre-buckling deformations are ignored. 

(f) The edges of each element along x = 0 and x = a, 
figures 23 and 26, are simply supported in the 
classical sense. 

The analysis has been coded for the CDC 6600 computer. 

The computer program "BUCLASPS", (Ref. 10) permits linear or 
parabolic thermal gradients across the width of individual 
plate-strip elements. For buckling analysis, the program has 
an input option whereby the user may also input desired thermal 
loads as preloads in the structure. The results quoted in this 
report are from the program "BUCLASP3." 
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4.0 THERMAL STRESS RESULTS 


The numerical results from the thermal stress analysis 
of Appendix B are correlated here with other analytical results. 
The choice of examples for correlation is restricted by the 
lack of published results for laminated composite structures 
including stiffened panels. However, the limited results 
presented here serve the basic purpose of verifvinq this analysis 
and do not reflect its full capability. 


The analysis assumes for each element simply supported 
boundary conditions along the edges normal to the longitudinal 
axis of the structure. The idealization of the structure as 
in the present analysis appropriately takes into account the 
relative internal stiffnesses and restraints which obviously 
influence the thermal stress distribution. All numerical results 
from the present analysis quoted below were obtained using the 
associated computer program "BUCLASP2", Ref. 10. 

(i) Isotropic flat panel with two "rod" stiffeners . -Figure 

1 shows the geometry of a rectangular isotropic flat panel with 
two stiffeners. The stiffeners are treated as rods (i.e., they 
have extensional stiffness only) and are symmetric with respect 
to the plate. A parabolic thermal gradient corresponding to 
values of a T AT of 0.0001 along the sides of the panel and 

0.0037 along the line of symmetry, is assumed. The rods are 
at a constant <* T AT of 0.0033. The results of the present 

analysis retaining 16 terms in equation (B.2) are shown in 
Figure 2 and 3. Plate numbers (1) and (2) refer to the 
corresponding numbers shown in Figure 1. For comparison the 
results from a finite element program. Ref. 11, are also shown. 

(ii) Isotropic cylinders . -The three isotropic cylinders of 
Figure 4 have constant temperature rise in the axial direction. 
The circumferential variation of temperature rise in each case 
is as indicated in the figure. 


The results from the present analysis are shown in Figures 
5 to 7. These results were obtained retaining 26 terms in 
equation (B.2). In Figure 5, the results for cylinder A are 
shown superposed on the results from Ref, 12. Some discrepancy 
in the N 22 value near the ends of the cylinder (x = 0) is 
noticed. This arises from the fact that in the present analysis, 
for each Fourier harmonic m in equation (B.2) the corresponding 
N 22 ^ as a " s i ne " distribution in the axial direction. Thus, 
as the number of harmonics is increased the results of the 
present analysis in the limit will approach the results of Ref. 
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For cylinder B, the thermal stress-resultants are shown 
in Figure 6 superposed on the results from Ref. 13. For both 
the cylinders A and B the thermal displacements, though not 
given here, showed excellent correlation. 

The results for cylinder C are shown in Figure 7 superposed 
on the results of Ref. 5. 

(iii) Laminated cylinder . -Figure 8 shows a laminated cylinder, 
each lamina being of a different isotropic material. The 
laminate has non-zero elements in the coupling matrix B of 
equations (A. 5) and (A. 6) . The cylinder is subjected to a 
uniform temperature rise of 400° F. The thermal stress-resultants 
are shown in Figure 9. The results from the present analysis 
are obtained retaining 16 terms in equation (B.2). Also shown 
are the results from Ref. 13. The cause of the discrepancy 
between the two results near the end of the cylinder (x = 0) 

is the same as that for the isotropic cylinder A. A3 the number 
of terms retained in equation (B.2) is increased the results 
of the present analysis will converge towards those of Ref. 

13. 

(iv) Stiffened cylinder- An integrally stiffened cylinder is shown 
in Figure 10. The stiffeners are idealized as beam elements 

and the parts of the shell between stiffeners as curved 
plate-strip elements. For a uniform temperature rise of 
AT = 100° F, the results of the present thermal analysis 
retaining 26 terms in equation (B.2), are shown in Figures 
11 to 13. Also shown are the results obtained from the 
computer program "BAMSOC 1" of Ref. 4, using classical analysis. 

To facilitate comparisons, the results quoted from the present 
analysis were obtained suppressing the underlined terms in 
equations (A. 4) and (A. 12) to (A. 14). The analysis of Ref. 4 
while taking into account the stiffener eccentricity, "smears" 
the properties of the stiffener. A similar method is used 
in Ref. 5. The present analysis not only retains the dis- 
creteness of each stiffener but also retains the relative 
stiffener (especially for other types of stiffeners) using 
the idealization discussed in Appendix A. The integral 
stiffeners of the cylinder under discussion were idealized as 
beam elements purely to conform with Ref. 4. It is more 
appropriate to idealize these stiffeners as flat plate-strip 
elements. 

The result of the present analysis reflect the effect 
of the discreteness of the stiffeners. Figure 11 shows the 
expected change in the out-of-plane displacement w of the 
shell skin in the region of each stiffener. Similar effects 
are noticed in the thermal stress distribution in the shell 
skin shown in Figure 12. The hoop stress N 22 from the 

present analysis is seen to be compressive everywhere in the 
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cylinder, whereas the results from " BAMSOC 1" show tensile 
hoop stress in the middle region of the cylinder. This 
difference is thought to be due to the discreteness of the 
stiffeners. For the equivalent "smeared" cylinder the present 
analysis will yield hoop stress distribution as in Figure 5 
for the unstiffened cylinder A, this being similar in nature 
to "BAMSOC 1" results. The peak in H 22 from the present analysis 

will move closer and closer to the cylinder end (x = 0) together 
with a simultaneous decrease in the unevenness of N 2 2 in the 

axial direction between x = 0.5 and x = 1.5, as the number 
of terms retained in equation (B.2) is increased. Another 
interesting feature of the results is the presence of the 
axial stress-resultant Njj . For a simply supported cylinder 
as in this example, analyses based on "smearing" the stiffeners 
do not predict this stress-resultant. Figure 13 shows the 
axial load find the bending moment in the integral stiffener. 

Similar results may be obtained for cylinders with different 
types of stiffening and also when circumferential thermal 
gradients are present. 
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5.0 BUCKLING RESULTS 


Numerical results of the buckling analysis of Appendix C 
are correlated here with other analytical results. 

As for the thermal stress results, the choice of examples for 
correlation is restricted by the lack of published results for 
heated laminated composite structures. The results presented 
here are intended to verify the analysis and does not indicate 
its full capability. 

(i) Rectangular flat plates . -Figure 14 shows a series of flat 
plates under different loading conditions. Each plate is 
idealized as an assemblage of four plate-strip elements. In 
plate A and B the only thermal stress is inplane and normal 
to the x axis. This stress, N 22 is uniform in the plates, 
thus making it easy to verify the results from the present 
analysis. 

From Ref. 14, the critical value of N 22 for plate A is 
1530 lbs/in at m = 1. The corresponding temperature rise 

N 2 2 

(AT =-p — - — — ) is 153° F. The present analysis yields at 
11 T z 

at buckling, a temperature rise of 155° F and n = 1. 

Plate B is subjected to uniform temperature rise of 35° F 
and axial mechanical loading Njj . The value of N 2 2 for 

this temperature rise is 350 lbs/in. The critical value of 

N 1;l from the present analysis is 2338.2 lb/in compared 

to 2334.8 lbs/in_from Ref. 14 for a biaxially loaded plate 
with a constant N 22 of 350 lbs/in. 

Plate C is subjected to uniform inplane bending as shown 
in Figure 14. Each plate-strip element is further divided 
into sub-elements. For a unit bending moment the linearly 
varying bending stresses are approximated by the corres- 
ponding averaged value in each sub-element. These values 
are increased by a constant factor until the plate buckles. 

From the present analysis the maximum bending stress-resultant 
N n at buckling is 272.5 lbs/in at m = 2 compared to 271 lbs/in 

and the same m value from Ref. 15, 

Plate D is subjected to a constant inplane bending 
corresponding to a maximum stress-resultant N 0 ** 46.885 lb/in 



with a superposed axial compression of II n as in Figure 14 

The bending stresses are considered as_preloads due to thermal 
stresses and the axial compression of II n is considered as 

further mechanical loading. The at buckling is obtained 

as 47.3 lbs/in at m = 2 compared to 46.9 lbs/in, from Ref. 15, 
at the same m value. 

(ii) Isotropic cylinders . -The present buckling analysis assumes 

the thermal stresses to be invariant in the axial direction. 
However, as stated in Appendix C, stresses averaged over 
the length nay be used when the axial variation is not too 
severe. Figure 5 shows the hoop stress resultant N 22 for 

cylinder A of Figure 4. It is seen that N 22 has a severe 
variation in the longitudinal (x) direction. Thus, the 
results from the present analysis can be expected to be 
very conservative. For this cylinder, Ref. 12 shows that at 
buckling the maximum value of the compressive hoop stress 
resultant No? near the end is about 3.44 times the value of 

the axially invariant hoop stress resultant causing buckling. 
Based on the true axial variation of U 2 2 • Ref. 12 shows the 

temperature rise at buckling for cylinder A to be 2380° F. 

The present analysis using averaged II? 2 values yields the 

critical temperature rise as 740° F. The ratio between these 
temperatures is 3.22, compared to 3.44 from Ref. 12. 

Cylinder C of Figure 4, heated on an axial strip is 
considered next. The cylinder is idealized by a aeries of 
curved plate-strip elements, the elements being considerably 
narrower in the region of heating. The thermal stress 
distribution in this cylinder for AT 0 = 1°F is shown in 
Figure 7. For this cylinder, the present buckling analysis 
gives the critical temperature { AT 0 ) = 286° F. The 

buckling mode shape is shown in Figure 15. This mode 
shape agrees with Ref. 5, where as the critical temperature 
quoted there is ( AT Q ) cr = 194° F. The buckling results of 

Ref. 5 is based on treating the cylinder as a portion of 
a very slender torus. This modelling is equivalent to 
considering an infinite cylinder, resulting in the thermal 
stress resultant Np to be constant in the axial direction. 

In the circumferential direction the value of H-q is close 

to the "two ring" case shown in Figure 6 of Ref. 5. Thus, the 
buckling analysis of Ref. 5 in effect approximates the axial 
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variation of Hu in the finite cylinder by its maximum value 
in the longitudinal direction compared to the averaged 
value used in the present analysis. Doth buckling analyses take 
into account the circumferential variation of . From Ref. 5 
the maximum thermal stress resultant corresponding to the 

temperature rise (AT 0 ) = 1°F, 

is 1.3 lbs/in compared to the average value of 0.89 lbs/in 
from the present analysis. Hence it is reasonable to expect 
the critical temperatures from the two analyses to be in the 
inverse ratio of these stress resultants. This in fact is 
seen to be the case. 

The cylinder D shown in Figure 1C is investigated 
for buckling under pure bending. This is equivalen to a 
case where the thermal stresses vary in the circumferential 
direction only. Analytical solutions for such cylinders 
are given in Refs. 16 and 17. 

In the present analysis cylinder D is idealized as a 
series of curved plate-strip elements which are further divided 
into sub-elements. These sub-elements being of sufficiently 
small width, the bending stress in each of them is approximated 
by an equivalent average stress. These average stresses in each 
sub-element corresponding to a unit bending moment are input 
into the computer program "BUCLASP3" , Ref. 10. The bending 
moment is increased in steps proportionately scaling the 
bending stresses until the cylinder buckles. This is equiva- 
lent to a thermal buckling problem. The results including the 
mode shape are shown in Figure 17. The buckling loads are 
given as the maximum compressive stress resultant Mji corres- 
ponding to the bending moment at buckling. Two sets of results 
from the present analysis are given. One set corresponds to 
the complete equations (A. 4) and (A. 12) to (A. 14) while the other 
set stems from reducting these equations to Donnell-type by 
suppressing the underlined terms. For comparison, the results 
from a computer program based on Ref. 18 are also given. The 
latter is based on Donnell-type equations. 

Figure 18 is a plot of the results from the present 
analysis for the buckling of cylinder D under pure bending to- 
gether with the results for buckling under uniform axial 
compression. The range of m-values studied is 1 to 19. 

As the m-value is increased the buckling load under bending is 
seen to approach the buckling load under axial compression. 

These results confirm the conclusions of Refs. 16 and 17. 

Similar results may be obtained for cylinders under 
arbitrarily varying axial and hoop stresses, the variation 
being confined to the circumferential direction. The 
cylinders may be of artitrary uniform cross section and may also 



be stiffened in the axial direction. 

The idealization of the structure as discussed in 
Appendix A retains the exact discreteness of the stiffeners 
and enables to determine the lowest of all buckling loads 
irrespective of whether buckling is of a general nature or is 
confined to a local region in the stiffener (or shell) . The 
eigenvector solution identifies the true nature of buckling. 

(iii) Stiffened cylinder .-Figures 12 and 13 show the thermal 
stress distribution in the integrally stiffened cylinder of 
Figure 10. These stresses exhibit severe variation in the 
axial direction. A3 discussed earlier, in such cases the 
present buckling analysis based on the thermal stresses averaged 
in the axial direction could yield very conservative solutions. 
Ignoring this for the moment, the stiffened cylinder under 
consideration is treated as an example where the thermal stresses 
do not vary in the axial direction, while the variations in the 
other directions are retained. 

The critical temperature rise at buckling from the present 
analysis is 238.4° F. The same temperature, on suppressing the 
underlined terms in equations (A. 4) and (A. 12) to (A. 14) , is 
247.7 F. In either case buckling occurs at n = 1. Thus, for 
this cylinder the effects of the additional terms are not very 
significant, for this particular buckling mode. The node shape 
is shown in Figure 19. The stiffeners are seen to enforce 
buckling nodes in the shell skin. Obviously, the "smearing" 
of the stiffeners as in Refs. 4 and 5 cannot yield such a 
buckling mode. For reasons stated in Section 4.0, the stiffeners 
in the cylinder have been idealized as beam elements. Such 
an idealization is not capable of yielding any possible local 
buckling in the stiffeners. This drawback is easily overcome 
by idealizing the integral stiffeners as flat plate-strip 
elements. The results for the cylinder under discussion thus 
brings out the capability of the present analysis in not only 
retaining the discreteness of the stiffners, but also their 
true structural characteristics. 

(iv) Sheet-stiffener unit . -A titanium sheet-stiffener unit 
typical of a supersonic airplane wing panel is shown in 
Figure 20. Here, the thermal gradients corresponding to 
descent conditions are approximated by constant temperatures 
in each segment. The axial thermal stress-resultants Njj in 

in these segments, calculated on the basis of uniform axial 
strain in the cross section are tabulated in the figure. 

The results of two simplified studies illustrating 
another application of the present analysis to buckling under 
axial compression load, are given in Figure 21. The riveted 
connection between the stiffener and the sheet is idealized as 
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in Refs. 8 and 10. The material properties used in corres- 
ponding to the temperature in each segment are given in 
Figure 20. In the first study, the sides A and B have 
boundary conditions corresponding to symmetry. The minimum 
buckling load occurs at m = 1. The resulting mode shape 
is interesting and shows coupling between Euler, torsional 
and local modes of buckling. In the second study the sides 
A and B are simply supported. The minimum buckling load in 
this case occurs at m = 18 and the corresponding mode shape 
indicates local buckling in the stiffener web. 

It may not be possible to predict the behavior of an 
entire stiffened panel from such results for an isolated sheet- 
stiffener unit. This is confirmed by the buckling mode shape 
for a zee stiffened panel given in Ref. 7. 



6.0 CONCLUDING RDI1ARKS 


Linear analyses for thermal stresses and buckling of 
structures of arbitrary uniform cross section have been 
presented. The structure is idealized as an assemblage of 
laminated curved plate-strip elements, laminated flat plate- 
strip elements and beam elements, each element extending the 
entire length of the structure. The element edges normal to 
the longitudinal axis are assumed to be simply supported with 
no restriction on the axial (warping) displacement. The 
idealization permits differentiating between bonded and riveted 
connections. The analyses consider the effects of offsets 
betx^een elements. Longitudinally stiffened panels and shells 
are typical examples of such structures. 

The idealization is seen to represent the true 
characteristics of these structure in a detailed manner. The 
relative internal stiffnesses and restraints which obviously 
influence the thermal stress distribution are appropriately 
represented unlike the more common method of either smearing 
the stiffeners or treating them as beams along their elastic 
axis. 


The buckling analysis ignores all pre-buckling deformations 
and treats the inplane normal thermal stresses as preloads. 

The preloading is assumed to be constant in each element in 
the longitudinal direction, while varying across the width of 
each plate-strip element. Thus the buckling analysis is also 
applicable to stiffened structures under inplane loads varying 
through the cross section as in stiffened shells under bending. 
The theory is general and no assumption is made regarding the 
buckling mode. The structure is free to take the buckled shape 
corresponding to minimum energy conditions consistent with 
prescribed constraints along any external side of plate-strip 
elements. The eigenvector solution is used to determine the 
buckled mode shape. The resulting mode shape plot is useful 
in understanding the buckling mechanism and may be used in 
efficient design of structures prone to buckling. 

The analyses assume that the material is linearly elastic. 
However, the elastic properties corresponding to the average 
temperature in each element may be used, 

A general lack of other published results makes it difficult 
to correlate the results from the present analysis. The limited 
results presented show good correlation. 
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APPENDIX A 


BASIC CONSIDERATIONS 


This appendix outlines the basic considerations involved 
in the analyses of Appendices B and C for thermal stresses and 
buckling, respectively. A summary of the basic equations used in 
the analyses are given. These equations are the result of in- 
cluding thermal effects in the basic equations of Ref. 7. 


A.l Idealization of the Structure 

Any structure of uniform arbitrary cross section is 
idealized as an assembly of laminated, curved and planar, plate- 
strip elements, and beam elements. The beam elements are used 
to idealize lips and beads in structural sections or any local 
reinforcement in the form of a lumped area of material. Each 
element extends the full length of the structure. The 
intersecting angle between elements is arbitrary. Element edges 
normal to the longitudinal axis are assumed to be simply 
supported with no restriction on the axial (warping) 
displacement. Arbitrary elastic restraint boundary conditions 
may be specified along any external longitudinal side (not 
connected to other elements) of flat and curved plate-strip 
elements. 

Figure 22 shows the uniform cross section of an arbitrary 
structure. The longitudinal (x) sides along G and II are 
arbitrarily restrained. This structure is idealized as an 
assemblage of: 

(a) Curved plate-strip elements: 2, 3 and 9 

(b) Flat plate-strip elements: 1, 4, 5, 7, 8, and 11 

(c) Beam elements : 6 and 10 

In this figure the element numbers are shown circled. Each 
plate-strip element is considered replaced by a membrane in 
its mid-plane (chosen as the reference plane) . Each beam element 
is considered replaced by a line through its shear center. 

The broken line shown is drawn through the mid-plane of each 
plate-strip element. The extremities of the individual mid- 
planes and the shear center of each beam element are marked 
with "dots." x Q , y and z G are 

the chosen global axes; x G being parallel to the longitudinal 

axis of the structure, y and z are the local axes in each 
element. The local x axis (not shown in the figure) is 
parallel to x Q . The idealization of the 



structure, as indicated above, also defines any offsets between 
elements. It is possible to differentiate between bonded and 
riveted connections as discussed in Ref. 10. 


A. 2 Laminated Curved and Flat Plates 

The equations given here are for the laminated curved 
plate. They degenerate to those of the laminated flat plate 
when the curvature becomes zero (infinite radius). 

Figure 23 shows the geometry and sign conventions for 
the curved laminate. The x, y and z axes are assumed 
coincident with the fiber axes 1, 2, and 3. The mid-plane of 
the laminate is chosen as the reference plane. The stress- 
strain equations for an orthotropic lamina with a temperature 
change T from ambient conditions (in the unstressed state 
T = 0) are, Ref. 19: 
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a Tx a Ty are t * ie coe ff^ c *- ents °£ thermal expansion of the 

layer in the two orthogonal direction x and y, respectively. 

Based on the Kir chhoff -Love hypothesis , the strains in 
any plane at distance z from the reference plane are: 




The above equations reduce to the Donnell-type assumptions when 
the underlined terms are dropped. 
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Combining equations (A. 1) , (A. 3) , and (A. 4) and integrating 
over the thickness of the laminate, the stre33 resultants and 
the moment resultants in the reference plane are: 
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h^ and h^ + ^ are the distances to the upper and lower surfaces, 

respectively, of the lamina. The A, B and D coefficient 
matrices define the overall extensional, coupling and bending 
stiffnesses, respectively, of the laminate, in relation to the 
chosen reference plane (here, the mid-plane). The thermal terms 
are shown with an additional subscript T. 
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Suitable "stability equations" and consistent boundary 
conditions for the_laminated curved plate under external biaxial 
mechanical loads, and M 2 2 / are derived in Ref. 7, ignoring 

pre-buckling deformations. The equations are: 
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After substituting equations (A.4) to (A. 6), the above equations 
are written as: 
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where the linear differential operators L i to L, are defined 
by: 5 
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Similarly, along any edge x — constant; 
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Q, Q, etc., are the effective force resultants along the side3 
and edges of the plate-strip element. The above equations 
(A. 12) to (A. 31) reduce to Donnell-type equations for thin 
cylindrical shells when the underlined terms are dropped. 
Further, when the external mechanical loads II ^ and N 2 2 

are zero, the above equations are identical to the equilibrium 
equations and the boundary conditions for thermal stress 
analysis. 

The displacements and forces in equations (A. 24) to (A. 31) 
are mid-plane values (chosen as the reference plane) and are 
with respect to the local axes. The positive directions of 

such displacements and forces along the sides y = + - 

and y = - - of the plate-strip element are shown in Figure 

24. Later, in Appendices B and C it will be necessary to 
transform these to parallel offset planes and also to global 
axes which make an arbitrary angle with the local axes. These 
transformations are discussed in Ref. 7 and are outlined 
below. 


Figure 25 shows a plate-strip element ABCD in relation 
to arbitrarily chosen global axes Xq, y^ and z Q . Also shown 
are typical offsets to 5, for the side nc. Axe3 system x g , y g 

and z g at S are chosen parallel to the local axes at B. The off- 
sets are defined by the distances y D and z Q measured positive 

along the positive local y and z axes, respectively, from 
the side of the plate-strip element in the reference plane. 

The displacements and forces in equations (A. 24) to (A. 27) are 
transformed to the axes system through S as: 
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Figure 25 also shows positive directions of the global 
displacements and forces, subscripted G. From geometry, it is 
seen that the transformation angle to global axes at A is: 

<6=^-1 (A. 40) 

A 2 

and at B 

0 B = t + \ (A. 41) 

The angle is measured positive in the clockwise direction, 
from the global axes. For a flat plate-strip element 


A matrix T Q for coordinate transformation through angle is 


defined as: 
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force vectors as: 
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and using superscripts - and + to denote the two sides of 
the plate-strip element identified by the values of 

b ^ b 

y - -- and y - +- , respectively, the displacements and 



forces in equations (A. 32) to (A. 39) are transformed to global 
axes as : 



(A. 47) 


(A« 4 8) 


T a and T D are the coordinate transformation matrices obtained 

from equation (a. 42) for angles of equation (A. 40) and 

of equation (A. 41) , respectively. The negative sign associated 
with M g in equation (A. 45) facilitates the use of a common 

coordination transformation matrix. The negative sign associated 
with T a in equation (A. 48) is a consequence of the sign con- 
vention for forces in the local axis system. 


A. 3 Laminated Beams 


Concentrated local reinforcements, like beads or lips in 
stiffeners, corner fillets in extruded sections, any beam type 
boron reinforcements, etc., are idealized as beam elements. 
Figure 26 shows the geometry and sign conventions for these 
elements. Any beam element is idealized as a line in the 
longitudinal (x) direction through its shear center. 


The general theory, Ref. l,_for beams subjected to a tempera- 
ture change T and axial load P fa , gives the following equations: 
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The second term on the right-hand side of equation (A. 56) 
represents the thermal effect on torsional stiffness. The presence 
of the term can either decrease or increase the effective 
torsional stiffness, though generally a decrease will occur. 

This is discussed in Refs. 1, 21, and 22. 

The positive directions of the displacements u, v, and 
w and the twist 6 of the beam element are shown in' Figure 26 . 
y^ and z m are the distances measured parallel to the 
principal axes y and z, respectively, from the shear center 0 
to the neutral axis C. These distances are positive in the 
positive directions of the axes. Approximate methods to evaluate 
the gross beam properties Ejjlyy, E nl zz » E n alp , C, 2 yJ 

and EjjAjj for the more common beam types, namely the laminated 
rectangular and the laminated circular beam elements are given 
in Ref. 7. 

Some limitations of the above equations for the heated beam 
element are discussed in Ref. 23. However, the equations given 
are considered to be sufficiently accurate for the envisaged 
applications of the beam element, mentioned at the beginning 
of this sub-section. 

An offset transformation identical to that for the 
plate-strip elements, equations (A. 32) to (A. 39) can be 
readily done. However, it is not pursued here as the offsets 
in plate-strip elements would suffice for the idealization 
of the types of structures considered. The beam element 
displacements and corresponding forces may be transformed 
to global axes oriented at an angle 0 b measured clockwise 

from the global axes to the local axes. Thus, 
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^bG anc * ^bG are global displacement and force vectors, 
defined in equations (A. 44) and (A. 46). T b is the coordinate 
transformation matrix of equation (A.42) , for the angle 4> h . 
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APPENDIX B 


THERMAL STRESS ANALYSIS 


The equations of Appendix A are now used in developing an 
analysis for thermal stresses in structures of uniform cross 
section. The structure is idealized as discussed in sub-section 

A. 1. 


For each plate-strip elements, linear or parabolic thermal 
gradients may be specified across the width. An arbitrary 
thermal gradient through the thickness is approximated by a 
series of layers, each at constant temperature. Temperature 
is constant in the longitudinal (x) direction. For the laminated 
circular and rectangular beam elements the temperature is 
arbitrarily constant in each lamina and in the longitudinal 
direction. For the latter type of beam element, identically 
linear or parabolic thermal gradients are allotted across the 
width of each lamina. The equations of sub-section A. 3 are 
applicable to any arbitrary beam element, provided the necessary 
beam properties are known. 

The material of the structure is assumed to be linearly 
elastic. However, the elastic properties corresponding to the 
average temperature in each lamina may be used. 

The constant temperatures in the axial direction of each 
element are approximated by a truncated Fourier sine series. 

The governing differential equations for each type of element 
are satisfied in an exact manner. Then, the stiffness matrices 
for these elements, relating the forces and displacements along 
inter-element boundaries, are derived taking into account the 
effects of thermal gradients. 

For each harmonic of the Fourier series, the element 
stiffness matrices are merged in a manner analogous to the 
direct stiffness method of finite element analysis. The axial 
(x) variations of the various quantities are separated out. 

The resulting non-homogeneous equations are solved for the 
thermal displacements using the symmetric Gaussian elimination 
algorithm. The displacements yield the thermal strains and 
stresses. Solutions are obtained in a similar manner for each 
haromonic and the results are superposed to yield the resultant 
thermal stresses. 


5.1 Stiffness of Heated Curved and Flat Plate-Strip Elements 

In a laminated curved or flat plate-strip elements, the 
temperature change from ambient conditions are specified as: 


T = T x 



(B.l) 
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Since the thermal field in the axial (x) direction, Figure 23, 
is assumed not to vary, T x is a constant for each plate-strip 

element and is approximated by the truncated Fourier series: 
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In the transverse (y) direction the thermal field is assumed as: 
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In the thickness direction (z) , the thermal gradient may be 
arbitrary and is approximated by a series of layers each at 
a constant temperature. Thus the temperature distribution in 

the k* 1 * 1 layer is given by: 

k k 

T = T • T • T 
x y . Z 


Equations are now developed for the stiffness matrix 
relating the thermal displacements and the corresponding forces 

i) 

along the sides y = ±- of the laminated curved plate-strip 

element. In the limit as the curvature becomes zero (infinite 
radius) , these equations degenerate to those of the flat plate- 
strip element. The displacements considered are: 


w, 



v and u 
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The corresponding forces are, from equations (A. 2 4) to (A. 27) : 
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Since the thermal stresses are being considered, all terms 
involving the external mechanical loads H and N-., have been 

dropped. The positive directions of the above displacements 
and forces are shown in Figure 24. The governing differential 
equations (A. 23) to (A. 25) for the plate-strip element are 
non-homogeneous . For each harmonic, m, in the equation (B.2) , 
the complementary solution is assumed as: 
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u = ^ ' U. e • 1 cos a 

c l—i 1 

i=l 


where a = “L and , = 

a H i a 

(i=l, 2, ...8) are the roots of the characteristic 

equation, discussed later. The corresponding particular solution 
is assumed as: 


/A A A 2\ . 

W P = n + V + K g y i sin “ 

/A A A 2\ 

V P = \ K 4 + K 5 Y + V ' Sin “ (B.7) 

/A A A 2 \ 

u p - + K 2 y + K^y / cos a 


The above displacement functions are chosen to satisfy ab initio 
the simply supported boundary conditions defined by: 
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along the edges x = 0 and x = a. For the complementary 
solution, a typical term of equation (B.6) is substituted into 
equations (A. 15) to (A. 17) and setting the right-hand sides to 
zero, results in: 



(B. 9) 


Expressions for the R^j coefficients (i, j *» 1, 2, 3) may be 

readily derived from equations (A. 18) to (A. 23) , after dropping 
all terms in Njj and H 22 • 0n expanding the determinant of 

the coefficient matrix in equation (B.9), a characteristic 
polynomial is obtained as: 

Vi 8 + Vi 6 + Vi 4 + Vi 2 + *o - 0 < B - 10) 

This equation yields eight values of p^ which are real or 

complex conjugates. Four of these roots are the negatives of 
the other four. Also, from equation (B.9) 
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The complementary displacement functions of equation (B * 6 ) 
are now fully defined for each harmonic, m, except for the 
amplitude coefficients W^. The contributions of these 

functions to the " complete” displacements and the corres- 
ponding "complete" forces of equations (B.4) and (B.5) may 
be readily calculated in terms of The procedure is 


analogous to that in Appendix B of Ref. 7. 

Next, the Kj_ coefficients (i = 1, 2, ...9) in the 

particular solution functions of equation (B.7) are deter- 
mined from equations (A. 15) to (A. 17) by the usual method 
for non -homogeneous differential equations. The right-hand 
sides in these equations contain all the thermal terms and 
are evaluated using equations (A. 10) and (A. 11^. The above 
procedure yields explicit expressions for the Kj_ coefficients. 

The details, though not given here, will result in all K £ 

A 

coefficients except Kg being non-zero. The particular 

solution functions are now fully defined. Thus, the 
contributions of the particular solution functions to the 
"complete" displacements and the corresponding "complete" 
forces of equations (B.4) and (B.5) are explicitly known. 


Finally, the complementary solution results, in terms 
of the yet unknown amplitude coefficients Wj , and the ex- 
plicitly known particular solution results are superimposed 
to yield the "complete" displacements and the "complete" 
forces for each harmonic, m, in equation (B.2). 


In idealizing a structure as discussed in Appendix A 
offsets between elements may occur. Further, the inter- 
secting angle between elements being arbitrary, it is con- 
venient to transform all displacements and forces to common 
global axes. The necessary transformations have been outlined 
in Appendix A. Following the procedure of Ref. 7, such trans- 
formations yield the following equations for the displace- 


ments and forces along the sides y 
the plate-strip element: 
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Offset transformation effects are included in the 
displacement and force vectors dp and fp, respectively, 

corresponding to the particular solution functions. dpQ 

and fpQ are the corresponding vectors after transformation to 

global axes. The matrices C^Xj and C 2 X 2 in conjunction with 

the vector give the displacements and forces, respectively, 

corresponding to the complementary solution functions. These 
matrices include offset effects. The matrices X 3 and X 4 are 

self-evident and result from transformation to global axes. 
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In the above equations all quantities except the vector 

are explicitly known. The matrix C 2 represents the axial 

(x) distribution of the displacements and forces. Since the 
Fourier harmonic m, is identical in all elements making up 
the structure, it is possible to separate out C 2 and drop 
it from further consideration. 

Substituting for Wj. from equation (B.14) in equation 

(B. 16) , the force-displacement relation for the laminated 
curved plate-strip element is written as: 
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All quantities in equations (B.19) except the displacement 
vector dg are explicitly known. The symmetric matrix [s] is 

analogous to the stiffness matrix of the plate-strip element 
derived in Ref. 7, Any arbitrary boundary conditions along 
the external longitudinal (x) side (not connected to other 
elements) of a plate-strip element are treated in the manner 
of Ref. 7. Thus, spring constants k w , k^, k v and k u with 
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reference to local axes x, y and z, specify the elastic 
restraints to freedoms corresponding to the displacements 
w , 0, v and u, respectively. A diagonal matrix k Q is defined 
as : 
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"Reduced" force-displacement equations may be derived from 
equation (B.19) by adopting the method of Ref. 7, when one 
side of the plate-strip element is elastically restrained as 

b 

above. Mien the side B (y = + - ) is elastically restrained 
the resulting equation is: 


f 


G 


s d 
A G 


f 

A 


where 


s 


A 


'11 




'21 
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+ f 
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f 

Ao 
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(B.23) 

(B.24) 

(B.25) 

(B.26) 

(B. 27) 


Similarly, when the side A (y 



is elastically restrained: 



(B. 28) 


where 
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22 


- s 


21 



'12 


(B. 29) 
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(B. 30) 


37 



+ 


(B . 31) 


f BO = - S 21 • ( S 22 + V) ' £ 


(B.32) 


The symmetric matrices s A and s 0 are analogous to the 

reduced stiffness matrices of the plate-strip element, derived 
in Ref. 7. 

All equations derived thus far are for the laminated curved 
plate-strip element. In the limit as the radius R -► °o , these 
equations degenerate to those of the laminated flat plate-strip 
element. For the latter case, some simplifications aiding 
numerical solution are possible, when the elements of the 
coupling matrix [b] in equations (A. 5) and (A. 6) are identi- 
cally zero. These are discussed in Appendix D of Ref. 7. 


B « 2 Stiffness of Heated Beam Elements 

The temperature change from the ambient conditions in 
beam elements has the same basic distribution as given by 
equation (B.l), As for the plate-strip elements, T x is 
assumed constant and is again approximated by the truncated 
Fourier series of equation (b.2). 

The displacement functions for a heated beam elements 

are : 


w = W, sin a 

D 

6 = ©, sin a 
b 


v = V, sin a 
b 


u = U, cos a 
b 


(B . 33) 


where a = — ^ — , m being the Fourier harmonic number. 

The above functions are chosen to satisfy the simply supported 
boundary conditions at the ends x = 0 and x = a. These 
displacement functions are substituted in the beam element 
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equations ^A.49) to (A. 52). For thermal stress analysis all 
terms in P^ and <7 are omitted since there are no external 
loads. This results in: 


f b = X 5 ^ + f b (B * 34) 

The vectors f t and d^ are defined in equations (A. 59) 

and (A. 60) , respectively. The elements of the matrix X 5 may 

be readily derived from equations (A. 49) to (A. 52). The force 
vector f^ is defined as: 


T 

(B. 35) 


and represents the thermal effects in the heated beam element. 

The elements of this vector are calculated for each Fourier 
harmonic , m, using equations (A.53) to (A. 55). As for the 
plate-strip elements, for each value of m it is possible 
to separate out and drop from further consideration, the axial 
distributions of the forces and displacements in equation (B.34) . 
Performing the transformation to global axes indicated in 
equations (A. 57) and (A. 58), the force-displacement equation 
for the heated beam element is written as: 


f b = 


d 2 M 


IT 


.2 

d M 


0 , 


L r=lv 


v 


ZT 
2 ' 


dx 


f bG * S b d bG + f bG 


where 


s = t • x • T, 
b b 5 b 


-1 


f = T f 
bG b b 


(B. 36) 


(B. 37) 


(B. 33) 


The symmetric matrix s^ is analogous to the beam element 
stiffness matrix derived in Ref. 7. 



B.3 Thermal Stress Analysis Formulation 

The equations developed thus far are now applied in deter- 
mining the thermal stresses for an arbitrary structure of uniform 
cross section shown in Figure 22. Sub-sections B.l and B.2 give 
the equations, for individual heated elements relating inter- 
element boundary global displacement vectors to corresponding 
force vectors for a chosen Fourier harmonic, m. These are: 

(i) "reduced" equations (B.23) and (B.28) for plate- 
strip elements with specified boundary conditions 
along one external longitudinal (x) side 

(ii) equation (B.19) for any other plate-strip elements 

(iii) equation (B.36) for beam elements. 

The above equations are given in matrix form. These 
equations are appropriately merged to satisfy inter-element 
compatibility and equilibrium in a manner identical to that 
in the direct stiffness method of finite element analysis. 

For each Fourier harmonic, m, the resulting equations may be 
written as : 


S D = F 


(B , 39) 


The symmetric stiffness matrix S is formed by appropriately 
merging the element stiffnesses, s A , s B and s^. Since there 

are no external mechanical loads, the vector F represents 
thermal effects and is formed from the vectors f r , f^, f B 

and f b . d is the yet unknown vector of inter-element boundary 

displacements corresponding to a particular m-value and is 
referenced to global axes. Equation (B.39) is solved for these 
displacements using the symmetric Gaussian elimination algorithm. 
The resulting d£ and djt vectors for each plate-strip 

element are substituted in equation (B.13) enabling the evaluation 
of the corresponding and hence the complementary solution 

functions of equation (B.6) . This, together with the already 
known particular solution functions of equation (B.7), defines 
the complete thermal displacement field for a specific value 
of the Fourier harmonic, m. The corresponding thermal strains 
in individual plate-strip elements are calculated from equations 
(A. 4). The thermal stresses in beam elements are easily cal- 
culated from equations (A. 49) to (A. 52). The thermal strains 
and the thermal stresses evaluated for each harmonic, m, 
are superposed to yield the resultant thermal strains and 
thermal stresses caused in the structure by the specified 
thermal field. The decreasing differences between successive 
resultant stresses, as the number of harmonics retained is 
increased, indicates that convergence is achieved. 
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APPENDIX C 


BUCKLING ANALYSIS 


In certain aerospace structural components like wing panels 
of high-speed aircraft, thermal stresses vary mainly in the 
chordwise and thickness directions while remaining nearly 
constant in the axial direction* Analogous situations exist 
in 3 ome shell type structures, Ref. 24, where the thermal 
stresses mainly vary in the circumferential and thickness 
directions. A linear buckling analysis of such heated structures 
of arbitrary uniform cross section is developed here. The 
structure is idealized as discussed in Appendix A. All 
prebuckling deformations are ignored. For buckling of these 
heated structures under mechanical loads, the inplane normal 
thermal stresses in plate-strip elements and the axial thermal 
load in beam elements are treated as preloads on the structure. 
The latter is kept constant while the former is increased until 
buckling occurs in some part of the structure. The thermal 
stresses for a given thermal gradient is directly proportional 
to the thermal level. Thus a buckling temperature for the 
structure under given thermal gradients may be determined by 
proportionately increasing the thermal stresses (for the given 
thermal gradients and for a chosen thermal level) , until buckling 
occurs in some part of the structure. 

The mechanical loading is basically axial compression. 
However, as in Ref. 7, selected plate-strip elements may have 
inplane biaxial mechanical loads, A typical example of the 
latter type of structure is heated stiffened panels under biaxial 
mechanical loads. In such panels it is reasonable to assume 
that the inplane mechanical loading transverse to the stiffener 
axis is carried entirely by the plate-strip elements in the 
plane of the skin. 

The buckling analysis of Ref. 7 is for structures where 
the loading within each plate-strip element is uniform. The 
method is versatile in that no assumptions are made regarding 
the buckling mode shape except that the half-wave length 

(-) of buckling is the same in all elements. Transversely, the 
ra 

structure is free to take the buckle shape corresponding to 
minimum energy conditions consistent with the constraints 
along any external longitudinal (x) side of plate-strip 
elements. The nature of buckling is determined from the plot 
of buckling mode shape resulting from the eigenvector solution. 

The method of Ref. 7 is extended here to allow for the 
variation of the preloads (thermal stresses) in the transverse 
direction of each plate-strip element. Because of these 



variations, the governing differential equations for each plate- 
strip element have variable coefficients. Using the method 
of Ref. 9, each such element is sub-divided and the variable 
coefficient differential equations are reduced to corresponding 
sets of equivalent constant-coefficient differential equations, 
one for each sub-division. 

Based on these reduced equations, a stiffness matrix for 
plate-strip elements relating the buckling displacement vectors 
along the inter-element boundaries to the corresponding forces 
is derived. For beam elements, the stiffness matrix derived 
in Ref. 7 is used. These element stiffness matrices are 
appropriately merged to yield the stiffness matrix of the total 
structure. The elements of this matrix are transcendental 
functions of the load level in each element and the longitudinal 

(x) half-wavelength (-) of buckling. For each value of m, 

m 

the buckling load is iteratively evaluated from the resulting 
non-linear eigenvalue problem in the manner of Ref. 7. The 
minimum of such buckling loads for all possible m-values 
is then the critical load of the structure. The corresponding 
eigenvector solution yields the mode shape. 


C.l Stiffness of Plate-Strip Elements with 
Transversely Varying Preloads 

Figure 27 shows the mid-plane (chosen here as the 
reference plane) of a curved plate-strip element, x, y and z 
are the local axes. The basic equations (A. 15) to (A. 17) for 
such elements include thermal effects. The terms with the 
additional subscript T are the resulting thermal terms. 

For the buckling analysis, these additional terms are to be 
dropped since the thermal stresses are treated as preloads. 
Subsequent references to these equations implicitly assume 
that these terms are dropped. The resulting equations are 
identical to those of Ref. 7 . These equations are based on 
uniform inplane normal external loads of Njj and 1 ^ 2 ' When 

the loads vary in the transverse (y) direction, it is easy to 
see that the resulting differential equations will have 
variable coefficients. This difficulty is overcome here by 
using the method of Ref. 9 in the manner outlined below. 

The inplane normal loads on the element are given by the 
stress resultants II ^ (y) in the x direction and II 22 (y) 

the y direction, both varying with y only. These are written 
as : 
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(C.l) 


Ni^y) 
^22 (Y) 


N lp (y) + N lm 


N 2 P (y) + N 2m 


(C. 2) 


Nip (y) and i^p (y) are the preloads (thermal stresses) 
varying with y only. N^ m and are the additional 
uniform mechanical loads. 


The plate-strip element is sub-divided into a series of 
equal sub-elements (1), (2), ....(g) as shown in Figure 27. 

As the number of sub-elements in increased, (thus decreasing 
their size) the amplitude of the variable coefficients in the 
differential equations decreases. In the limiting case, when 
the sub-elements are infinitesimally small, these variable 
coefficients approach a fixed value. The original differential 
equations with variable coefficients may then be approximated 
by sets of constant coefficient equations, identical to 
equations (A. 15) to (A. 17) , one set for each sub-element. A 
relatively small number of sub-elements will usually suffice 
to yield accurate results, as shown in Ref. 9. The constant 
coefficient equations for each sub-element are satisfied in 
an exact manner. Continuity is achieved between sub-element 
regions from compatibility and equilibrium considerations. 


Figure 7.1 also shows the sub-element (j) of width b and 
the local axes x, y and z. The assumed buckling displacement 
function for each sub-element is: 



a 


V 



h . 

e sm a 


i=l 


(C. 3) 


u 



i=l 



e cos a 


where : 

P. Try 

/J . = and cn 

l a 


m7rx 

a 


(i-1, 2, ..,.8) are the roots of the characteristic 
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equation discussed below. The longitudinal (x) half-wave 
number, m, is by necessity the same in all elements of 
the structure. The above functions are chosen to satisfy 
ab initio the simply supported boundary conditions along 
the edges x = 0 and x = a. At any level of external loading 
an average value of II (j) and N 22 (j) sub-element (j) 

may be calculated from equations (C.l) and (C.2). In practice, 
for each sub-element, the values calculated along y = 0 
would suffice. These external load values, together with 
buckling displacement functions, are substituted in equations 
(A. 15) to (A. 17) to yield a characteristic polynominal in p, , 
identical in form to equation (B.10). Using these 1 

values, the buckling displacements and the corresponding 

forces along the boundaries y = - ^ and y =» + - of each 

sub-element are evaluated as in Ref. 7. Thus, for sub- 
element ( j ) : 




C 1 X 2 


C 1 X 2 


J (j) 



(C. 4) 


(C. 5) 


These equations are with respect to the local coordinates. 
The elements of each sub-vector |d| and |f| are readily 
identified from equations (A. 24) to (A. 27). The matrices 
X and X are analogous to Xj and X 2 in equations (B.13) and 


(B.15) . Cj, defining the longitudinal (x) distribution of 


the forces and displacements, is given by equation (B.17). 
These distributions, being identical in all elements, may be 
separated out and omitted from further consideration. The 
above equations are used in ensuring continuity between the 
sub-elements. This may be achieved by two somewhat similar 
methods . 


In the first method a symmetric stiffness matrix is 
derived for each sub-element from equations (C.4) and (C.5) , 
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as in Ref. 7 , relating the forces and displacements. This is 
written as: 



(C.6) 


Now, considering sub-elements (1) and (2) , the respective 
sub-element stiffness matrices are appropriately merged and 
the degrees of freedom along the common boundary are reduced 
to give: 


( 1 ) 


( 2 ) 


/ - \ 

d (D 


= [s] 


( 1 - 2 ) 


(2) 


> 


(C. 7) 


where [s]^ ^ is t ^ ie symmetric reduced stiffness matrix. By 


successively merging the stiffness matrices of adjacent sub- 
elements and reducing the interior freedoms at each step, the 
resulting equation becomes: 


/_ 


f 1 

d 1 

(1) 1 

(D I 

r [s] (l-g) - 


_ A 

- + 1 

' f (g) ' 

d (g) ' 


where 


[ s l(l-g) 


is the required symmetric stiffness matrix of 


the curved plate-strip element, when the external loads are 
varying in the transverse (y) direction. Instead of suc- 
cessively merging and reducing at each sub-element level 
as described above, it is possible to merge the stiffness of 
all sub-elements at one time and perform one large reduction 
of all interior freedomes. 


In the second method, continuity conditions imposed through 
equations (C.4) and (C.5) are used to relate VI ^ s of successive 


>ub 


-elements to • Thus for sub-element 


w. 

ii 


( 2 ) 


or 
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-1 

“ — 
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X, 


X, 
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L 2 J 
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(i) 





(i) 


(2) t 


(C.9) 
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Similarly, 




( 2 ) 



(C.10) 


The matrices involved in the above equations are unsymmetric. 
Substituting for {w } from equation (C.9) : 

M2) 





(C.ll) 


The matrix [c] is self-evident. Continuing this procedure 
to the last sub-element (g) yields: 


w. 

A 


= [c] . , jw.l 
(g) (g) * ^ (l) 


(C. 12) 


Also, for the boundary y 


of sub-element (1) : 



(1) 


L x 2 J 


w. 

1 ' 


(C. 13) 


(D- 


( 1 ) 


Similarly for the boundary y = of sub-element (g) : 


or 



(C. 14) 


(C. 15) 


After partititioning suitably, the stiffness matrix of the 
curved plate-strip element under transversely (y) varying 
external loads is readily derived from equations (C.13) and 
(C.15). The result is identical to equation (C.8), 



In the first method the inversion process concerns symmetric 
matrices whereas in the second method unsymmetric matrices are 
involved. The stiffness matrix derived by either method is 
based on equations (C.4) and (C.5) which are in local coordinates. 
A coordinate transformation on equation (C.8) yields the global 
force-displacement relationship for a curved plate-strip 
element ass 



The sub-vectors d G and f G are defined in equations (4.44) 
and (A. 46) , respectively. 

The buckling analysis allots arbitrary elastic restraints 
(specified by spring constants) along any external longitudinal 
(x) side of the curved plate-strip elements. For such elements 
reduced stiffness matrices are easily derived as in Ref. 7. 

b 

Thus, when the side y = of the element is elastically re- 

strained: 



(C. 17) 


and when the side y =* - is elastically restrained: 


f + = s d + 
G B G 


(C. 18) 


s A and s B are the reduced stiffnesses. 

The elements of the stiffness matrices in equations (C.16) 
to (C.18) are transcendental functions of the load level in 
the structure and the longitudinal (x) half-wave number, m. 

The analysis presented in this section for curved plate- 
strip elements degenerates to that of a flat plate-strip element 
in the limit as the radius R -*■ «> . In this case certain 

simplifications aiding numerical solution are possible and 
are discussed in Appendix D of Ref. 7. 
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C.2 Stiffness of Beam Elements 

The axial thermal loads are treated as preloads. The 
external mechanical load is superposed on this. The stiffness 
of the beam element subjected to the resultant axial load is 
easily derived as in Ref. 7 and results in 


f 


bG 


s d 
b bG 


(C. 19) 


The vectors fv_ and d,-. are defined in equation (A. 53) 
and (A. 54). ^ ^ 


C.3 Buckling Formulation 

The equations developed thus far are applied in formulating 
the buckling criterion for the arbitrary structure of uniform 
cross section shown in Figure 22. The method being the same 
as in Ref. 7 , only an outline is given below. 

The structure is. idealized as described in Appendix A. 

By necessity the number of longitudinal (x) half-waves , 
m, is identical in all elements of the structure. As indicated 
earlier, the analysis assumes that the thermal stresses do not 
vary in the longitudinal direction. Ilov/ever, if such variations 
do exist and are not too severe, a buckling solution of 
reasonable accuracy may be obtained using thermal stresses 
averaged in the longitudinal direction. Alternatively, the 
peak stress in the longitudinal direction may be used as 
illustrated for the buckling of a cylinder heated on an axial 
strip, in Ref. 5. 

In buckling under combined thermal gradients and mechanical 
loads, the thermal stresses are treated as preloads which remain 
constant. Thus different levels of loading are defined by 
changes in mechanical loading. The distribution of mechanical 
loading among the various elements in the structure, 
corresponding to any level of loading, is obtained from strain 
compatibility considerations, as in Ref. 7. Superposing these 
element loadings on the preloads (thermal stresses) , the 
stiffness matrices are evaluated for various elements in the 
structure, for a chosen value of m. These are merged to yield 
the stiffness matrix of the total structure. This procedure 
is analogous to that described in Appendix B. The resulting 
equation is: 

s ' D = 0 (C. 20) 

where S is the merged stiffness matrix and n is the vector 
°f , the inter-element boundary displacements. For a non- 
trivial solution, the buckling criterion is v/ritten in the 
determinantal form: 

|s | = o (C. 21) 
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The elements of the symmetric matrix 5 are trans- 
cendental functions of the load level in the structure and the 

longitudinal (x) half-wavelength of buckling (-). Thus the 

m 

above equation does not correspond to standard algebraic 
eigenvalue problems. For the chosen longitudinal half-wave 
number, ra, the lowest level of external loading at which 
equation (C.21) is satisfied i3 the buckling load of the 
structure. This is determined iteratively using the algorithm 
described in Ref. 7. An upper bound to the buckling load of 
the structure, determined from the buckling load of individual 
elements when the inter-element boundaries are completely 
restrained, is crucial to the success of the algorithm. A 
method to evaluate the buckling loads of these "restrained" 
elements subjected to transversely (y) varying preloads and 
uniform mechanical loads, is discussed in sub-section C.4. 

The lowest of all buckling loads of the structure for a 
series of m-values is then the critical load. The corres- 
ponding vector D in equation (C.20) giving the inter-element 
boundary displacements, is obtained by Wielandt's method of 
inverse iteration. Ref. 25. 

In order to calculate the distribution of the buckling 
displacements across the width of each plate-strip element, 
it is necessary to know the coefficients in equation 

(C. 3) , for each sub-element. These are calculated from 
the known element boundary displacements, using equations 
(C.4) and (C.8) or equations (C.4) and (C.12) . 

For beam elements the eigenvector solution directly 
gives the global buckling displacements of the shear center. 

A plot of w and v displacements of each element 
across the cross section of the structure identifies the 
buckled elements and indicates whether buckling is local 
or general. Such plots may be used in achieving efficient 
design of structures as illustrated in Refs. 14 and 26. 

The generality of the buckling analysis presented here is 
discussed and illustrated in Ref. 7. 



C.4 Upper Bound to Buckling Load 

The principle of the method is to determine a buckling 
load P r for each element, corresponding to the chosen value 

of the longitudinal (x) half-wave number m for the structure, 
when the inter-element boundaries are completely restrained. 

The smallest of the P r values is then the upper bound load 

to the buckling load of the structure, for the chosen m-value. 
Obviously, the beam elements may be ignored from these con- 
siderations. Thus, only plate-strip elements with either edges 
fully restrained or one edge fully restrained and the other 
edge arbitrarily restrained by specified spring constants, 
are considered. 

In the manner of sub-section C.l, each "restrainted' 1 
plate-strip element is divided into sub-elements as shown in 
Figure 27. The transversely-varying preloads (thermal stresses) 
are approximated by uniform values in each sub-element. The 
buckling load of each "restrained" element (for a chosen m-value 
for the structure) is obtained by formulating the corresponding 
buckling problem in the form of equation (C.21) . This is 
achieved by merging the stiffness of each sub-element and 
taking into account the restraints along the longitudinal (x) 
sides of the "restrained" element. Since a completely re- 
strained side has zero displacements, the corresponding rows 
and columns of the merged stiffness matrix are deleted. The 
resulting matrix is the buckling determinant from which the 
buckling load is evaluated iteratively. An upper bound to the 
"restrained" element buckling problem is obtained using the 
Galerkin method, as described in Appendix E of Ref. 7. 
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Figure 2. Thermal Stress in the "Rod" (1 of 2) 
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Figure 6. Thermal Stress in Cylinder B (1 of 2) 
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Figure 7. Results for Cylinder C 
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Figure 8. Laminated Cylinder 
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Figure 9. Thermal Stress in Laminated Cylinders (1 of 2) 
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Internal Stringers 



Number of stringers = 30 

AT = 1 00° F 


Skin Properties : 

Ejj = E 22 = 10.5 x 10® lbs/in 2 
6^2 = 3. 98 x 10® IDs/in 2 

V]2 — ■ 0.32 

<Xj = 14.0 x 10"® i n/i n/°F 

Stringer Properties : 

GJ = 2.38 x 10 4 lbs/in 2 

I = 0.088 in 4 

E = 10.5 x 10® IDs/ i n 2 

a T = 14.0 x 10'® i n/i n/°F 


Figure 10. Stiffened Cylinder 
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Figure 15 Buckling Mode Shape for Cylinder C 

R = 29.948 inches 

t = 0.1039 inches 

Length = 94.25 inches 

Material Properties : 

E 11 = E 22 = 10 *° x 1()6 lbs / 1n2 
G 12 = 4.286 x 10 6 lbs/in 2 

^ 12 = 0.1667 

Figure 16 Geometry of Cylinder D 







Maximum Compressive Stress Resultant 

(lb/in) 


Present Analysis I* 

Present Analysis II** 

Ref. 18 

1 

2814.8 

2999.7 

2993.6 

2 

2609.9 

2700,2 

2690.1 

3 

2511.5 

2571.5 

2563.5 

4 

2453.7 

2499.1 

2489.5 

5 

2414.4 

2451.2 

2439.4 

6 

2390.8 

2422.1 

2402.3 

7 

2372.2 

2399.2 

2373.3 

8 

2351 .1 

2381 .0 

2349.7 

9 

2328.5 

2349.6 

2329.9 


* Retainning all terns in equations (4.4) and (4.12) to (4.14) 

** Dropping underlined terns in equations (4.4) and (4.12) to (4.14) 



Figure 17 Results for Buckling of Cylinder D Under Pure Bending 
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Present Analysis I 
Present Analysis II 


I 



Present Analysis II drops the underlined terns in the equations 

Figure 18. Buckling of Cylinder D Under Pure Bending and 
Under Axial Compression 



Figure 19. Mode Shape for Thermal Buckling of Stiffened 
Cylinder (m = 1) 
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Figure 24. Displacements and Forces Along the Sides of the 
Plate-Strip Element 



Figure 25. Offsets and Global Axes for Plate-Strip Element 





y = + b/2 


r / 



N ii(y) = N i P (y) + N im 

^22< y ) * ^2p(y) + ^2m 


Sub-element (j) 


Figure 27. Sub-elements in Curved Plate-Strip Element 




APPENDIX D 


CONVERSION OF U.S. CUSTOMARY UNITS TO SI UNITS 


The International System of Units (SI) was adopted by the Eleventh 
General Conference on Weights and Measures, Paris, October 1960, In Resolution 
No. 12. (See Ref. 27)) Conversion factors for the units used herein are 
given in the following tables: 


Physical quantity 

U.S. Customary 
Unit 

Conversion 

factor 

(*) 

SI Unit ! 

(**) | 

Area 

in 2 

6.452 x 10" 4 

2 1 

square meters (m ) 

Force 

kip = 1000 lbf 

4.448 x 10 3 

newtons (N) 

Length 

in. 

2.54 x 10“ 2 

meters (m) 

Moduli and stress 

ksl * 1000 lbf/ 
In 2 

6.895 x 10 6 

newtons per square 
meter (N/m 2 ) 

Stress resultant 

Ibf/in. 

175.1 

newtons per meter 
(N/m) 

Temperature change 

°F 

5/9 

Kelvin (K) 


* Multiply value given In U.S. Customary Unit by conversion 
factor to obtain equivalent value in SI Unit. 

** Prefixes to Indicate multiple of units are as follows: 


Prefix 

Multiple 

mllli (m) 

10“ 3 

centl (c) 

10’ 2 

kllo(k) 

10 3 

giga(G) 

10 9 
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